! include 'mkl_pardiso.f90'

module MySparseOperator
    implicit none
    integer, parameter :: mso_dp = 8    ! double presion



!     interface mso_sparse_transpose
!         module procedure mso_sparse_d_transpose, mso_sparse_z_transpose
!     end interface ! mso_sparse_transpose


contains


    subroutine mso_z_uptri(upper,a,ia,ja,b,ib,jb)
        ! return the upper triangular part of A
        ! refered to triu in matlab
        implicit none
        integer, intent(in) :: upper
        complex(kind=8), intent(in) :: a(:)
        integer, intent(in) :: ia(:), ja(:)
        complex(kind=8), intent(inout) :: b(:)
        integer, intent(inout) :: ib(:), jb(:)


        integer :: nnz, ii, jj, kk
        integer :: rows

        rows = size(ia) - 1


        ib = 0
        do ii = 1, rows, 1
            do jj = ia(ii), ia(ii+1)-1, 1
                if ( ja(jj)-ii.ge.upper ) then
                    ib(ii+1) = ib(ii+1) + 1
                end if
            end do
        end do
        ib(1) = 1
        do ii = 2, rows+1, 1
            ib(ii) = ib(ii-1) + ib(ii)
        end do

        nnz = ib(rows+1) - 1
        kk = 1
        do ii = 1, rows, 1
            do jj = ia(ii), ia(ii+1)-1, 1
                if ( ja(jj)-ii.ge.upper ) then
                    b(kk) = a(jj)
                    jb(kk) = ja(jj)
                    kk = kk + 1
                end if
            end do
        end do
    end subroutine mso_z_uptri



    subroutine mso_sparse_d_show(a,ia,ja, fmt, isfull)
        ! trans: 'N', 'n', 'T', 't', 'C', 'c'
        ! isfull: not pesent or >0 if full, <0 otherwise
        implicit none
        double precision, intent(in) :: a(:)
        integer, intent(in) :: ia(:), ja(:)
        character(len=*), intent(in) :: fmt
        integer, intent(in), optional :: isfull

        integer :: ii, jj, kk
        character(len=20) :: s1
        integer :: rows, cols, nnz

        rows = size(ia) - 1
        nnz = ia(rows+1) - 1
        cols = maxval(ja(1:nnz))

        if ( rows*cols.eq.0 ) return

        if ( .not.present(isfull) .or. isfull.gt.0 ) then
            do ii = 1, rows, 1
                if ( ia(ii+1)-ia(ii).eq.0 ) then
                    do jj = 1, cols, 1
                        write(*,'('//fmt//'$)') 0.d0
                    end do
                else
                    kk = ia(ii)
                    do jj = 1, cols, 1
                        if ( jj.eq.ja(kk) ) then
                            write(*,'('//fmt//'$)') a(kk)
                            if ( kk.lt.ia(ii+1)-1 ) kk = kk + 1
                        else
                            write(*,'('//fmt//'$)') 0.d0
                        end if
                    end do
                end if
                write(*,*)
            end do
            return
        else
            do ii = 1, rows, 1
                do jj = ia(ii), ia(ii+1)-1, 1
                    write(*,'(I5,I5,'//fmt//')') ii, ja(jj), a(jj)
                end do
            end do
        end if
    end subroutine mso_sparse_d_show



    subroutine mso_dcoo2csr(rowind,colind,acoo,acsr,ia,ja,m,nnz)
        integer, intent(in) :: m,nnz
        integer, intent(in) :: rowind(nnz), colind(nnz)
        double precision, intent(in) :: acoo(nnz)
        double precision, intent(inout) :: acsr(nnz)
        integer, intent(inout) :: ia(m+1), ja(nnz)

        integer :: job(8), info

        job(1) = 2  ! converted to CSR
        job(2) = 1  ! one-based indexing for CSR
        job(3) = 1  ! one-based indexing for COO

        job(6) = 3
        call mkl_dcsrcoo( job, m, acsr, ja, ia, nnz, acoo, rowind, colind, info)
        !         job(6) = 2
        !         call mkl_dcsrcoo( job, m, acsr, ja, ia, nnz, acoo, rowind, colind, info)
    end subroutine mso_dcoo2csr


    subroutine mso_zcoo2csr(rowind,colind,acoo,acsr,ia,ja,m,nnz)
        integer, intent(in) :: m,nnz
        integer, intent(in) :: rowind(nnz), colind(nnz)
        complex(kind=8), intent(in) :: acoo(nnz)
        complex(kind=8), intent(inout) :: acsr(nnz)
        integer, intent(inout) :: ia(m+1), ja(nnz)

        integer :: job(8), info

        job(1) = 2  ! converted to CSR
        job(2) = 1  ! one-based indexing for CSR
        job(3) = 1  ! one-based indexing for COO

        job(6) = 3
        call mkl_zcsrcoo( job, m, acsr, ja, ia, nnz, acoo, rowind, colind, info)
        !         job(6) = 2
        !         call mkl_dcsrcoo( job, m, acsr, ja, ia, nnz, acoo, rowind, colind, info)
    end subroutine mso_zcoo2csr



    ! old codes



    subroutine mso_zgmres_precondition(maxiter, maxrelres, a,ia,ja,rhs, computed_solution, itercount, relres, prefun)
        !*******************************************************************************
        !  Content:
        !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
        !*******************************************************************************
        implicit none
        interface
            subroutine prefun(x,u)
                complex(kind=8), dimension(:), intent(in) :: x
                complex(kind=8), dimension(:), intent(inout) :: u
            end subroutine prefun
        end interface
        integer, intent(in) :: maxiter(2)
        double precision, intent(in) :: maxrelres
        complex(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        complex(kind=8), dimension(:), intent(in) :: rhs
        complex(kind=8), dimension(:), intent(inout) :: computed_solution
        integer, intent(inout) :: itercount
        double precision, intent(inout) :: relres


        integer :: n, zn
        integer :: ipar(128)
        double precision :: dpar(128)
        double precision, dimension(:), allocatable :: tmp, b, drhs, dcs
        complex(kind=8), dimension(:), allocatable :: residual

        integer :: rci_request, i
        double precision :: dvar, rhsm2
        double precision, external :: dznrm2
        integer :: err

        err = 0

        zn = size(ia) - 1
        n = 2*zn    ! the size of the problem
        rhsm2 = dznrm2(zn,rhs,1)

        allocate( b(n))
        allocate( residual(zn))
        allocate( drhs(n))
        allocate( dcs(n))

        !---------------------------------------------------------------------------
        ! Initialize variables and the right hand side through matrix-vector product
        !---------------------------------------------------------------------------
        call dcopy(zn, dreal(rhs), 1, drhs, 1)
        call dcopy(zn, dimag(rhs), 1, drhs(zn+1), 1)

        !       dcs = 0.d0
        call prefun(rhs,computed_solution)
        call dcopy(zn, dreal(computed_solution), 1, dcs, 1)
        call dcopy(zn, dimag(computed_solution), 1, dcs(zn+1), 1)

        ipar(15) = maxiter(2)
        allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

        call dfgmres_init(n, dcs, drhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        ipar(8) = 1
        ipar(15) = maxiter(2)
        ipar(5) = maxiter(1)
        ipar(11) = 1    ! precondition
        dpar(1) = maxrelres

        call dfgmres_check(n, dcs, drhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        1 call dfgmres(n, dcs, drhs, rci_request, ipar, dpar, tmp)

        if (rci_request .eq.0) go to 3

        if (rci_request .eq. 1) then
            computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),    &
                    &   tmp(ipar(22)+zn:ipar(23)-1) )
        !           call zcopy(zn,dcmplx(tmp(ipar(22):ipar(22)+zn-1),tmp(ipar(22)+zn:ipar(23)-1)),1,computed_solution,1)
            call mkl_zcsrgemv('n',zn, a, ia, ja, computed_solution, residual)
            call dcopy(zn,dreal(residual),1,tmp(ipar(23)),1)
            call dcopy(zn,dimag(residual),1, tmp(ipar(23)+zn), 1)
            go to 1
        end if

        if (rci_request .eq. 2) then
            ! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
            ipar(13) = 1
            ! get the current fgmres solution in the vector b(n)
            call dfgmres_get(n, dcs, b, rci_request, ipar,dpar, tmp, itercount)
            ! compute the current true residual via mkl (sparse) blas routines
            computed_solution = dcmplx(b(1:zn),b(zn+1:n))
            call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
            call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
            dvar = dznrm2(zn, residual, 1)
            dvar = dvar/rhsm2
            if (dvar .lt. maxrelres) then
                go to 3
            else
                go to 1
            end if
        end if

        if (rci_request .eq. 3) then
            computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),    &
                    &   tmp(ipar(22)+zn:ipar(23)-1) )
            call prefun( computed_solution, residual )
            call dcopy(zn, dreal(residual), 1, tmp(ipar(23)), 1)
            call dcopy(zn, dimag(residual), 1, tmp(ipar(23)+zn), 1)
            go to 1
        end if

        if (rci_request .eq. 4) then
            if (dpar(7) .lt. 1.0d-12) then
                go to 3
            else
                go to 1
            end if
        else
            go to 999
        end if

        3 ipar(13) = 0
        call dfgmres_get(n, dcs, drhs, rci_request, ipar, dpar, tmp, itercount)
        computed_solution = dcmplx( dcs(1:zn), dcs(zn+1:n) )
        call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
        call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
        relres = dznrm2(zn, residual, 1)
        relres = relres/rhsm2

        if (allocated(dcs)) deallocate(dcs, stat=err)
        if (err /= 0) print *, "dcs: Deallocation request denied"
        if (allocated(drhs)) deallocate(drhs, stat=err)
        if (err /= 0) print *, "drhs: Deallocation request denied"
        if (allocated(tmp)) deallocate(tmp, stat=err)
        if (err /= 0) print *, "tmp: Deallocation request denied"
        if (allocated(b)) deallocate(b, stat=err)
        if (err /= 0) print *, "b: Deallocation request denied"
        if (allocated(residual)) deallocate(residual, stat=err)
        if (err /= 0) print *, "residual: Deallocation request denied"
        return


        999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
        call mkl_free_buffers
        stop 1
    end subroutine mso_zgmres_precondition

    subroutine mso_zgmres_nonprecondition(maxiter, maxrelres, a,ia,ja,rhs, computed_solution, itercount, relres)
        !*******************************************************************************
        !  Content:
        !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
        !*******************************************************************************
        implicit none
        integer, intent(in) :: maxiter(2)
        double precision, intent(in) :: maxrelres
        complex(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        complex(kind=8), dimension(:), intent(in) :: rhs
        complex(kind=8), dimension(:), intent(inout) :: computed_solution
        integer, intent(inout) :: itercount
        double precision, intent(inout) :: relres


        integer :: n, zn
        integer :: ipar(128)
        double precision :: dpar(128)
        double precision, dimension(:), allocatable :: tmp, b, drhs, dcs
        complex(kind=8), dimension(:), allocatable :: residual

        integer :: rci_request, i
        double precision :: dvar, rhsm2
        double precision, external :: dznrm2
        integer :: err

        err = 0

        zn = size(ia) - 1
        n = 2*zn    ! the size of the problem
        rhsm2 = dznrm2(zn,rhs,1)

        allocate( b(n))
        allocate( residual(zn))
        allocate( drhs(n))
        allocate( dcs(n))

        !---------------------------------------------------------------------------
        ! Initialize variables and the right hand side through matrix-vector product
        !---------------------------------------------------------------------------
        call dcopy(zn, dreal(rhs), 1, drhs, 1)
        call dcopy(zn, dimag(rhs), 1, drhs(zn+1), 1)

        dcs = 1.d0

        ipar(15) = maxiter(2)
        allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

        call dfgmres_init(n, dcs, drhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        ipar(8) = 1
        ipar(15) = maxiter(2)
        ipar(5) = maxiter(1)
        ipar(11) = 0
        dpar(1) = maxrelres

        call dfgmres_check(n, dcs, drhs, rci_request, ipar, dpar, tmp)
        if (rci_request .ne. 0) go to 999

        1 call dfgmres(n, dcs, drhs, rci_request, ipar, dpar, tmp)

        if (rci_request .eq.0) go to 3

        if (rci_request .eq. 1) then
            computed_solution = dcmplx( tmp(ipar(22):ipar(22)+zn-1),    &
                    &   tmp(ipar(22)+zn:ipar(23)-1) )
            call mkl_zcsrgemv('n',zn, a, ia, ja, computed_solution, residual)
            call dcopy(zn,dreal(residual),1,tmp(ipar(23)),1)
            call dcopy(zn,dimag(residual),1, tmp(ipar(23)+zn), 1)
            go to 1
        end if

        if (rci_request .eq. 2) then
            ! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
            ipar(13) = 1
            ! get the current fgmres solution in the vector b(n)
            call dfgmres_get(n, dcs, b, rci_request, ipar,dpar, tmp, itercount)
            ! compute the current true residual via mkl (sparse) blas routines
            computed_solution = dcmplx(b(1:zn),b(zn+1:n))
            call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
            call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
            dvar = dznrm2(zn, residual, 1)
            dvar = dvar/rhsm2
            if (dvar .lt. maxrelres) then
                go to 3
            else
                go to 1
            end if
        end if

        if (rci_request .eq. 3) then
            write(unit=*, fmt=*) 3

            call dcopy(n,tmp(ipar(22)),1, tmp(ipar(23)), 1)
            go to 1
        end if

        if (rci_request .eq. 4) then
            if (dpar(7) .lt. 1.0d-12) then
                go to 3
            else
                go to 1
            end if
        else
            go to 999
        end if

        3 ipar(13) = 0
        call dfgmres_get(n, dcs, drhs, rci_request, ipar, dpar, tmp, itercount)
        computed_solution = dcmplx( dcs(1:zn), dcs(zn+1:n) )
        call mkl_zcsrgemv('n', zn, a, ia, ja, computed_solution, residual)
        call zaxpy(zn, -1.0d0, rhs, 1, residual, 1)
        relres = dznrm2(zn, residual, 1)
        relres = relres/rhsm2

        if (allocated(dcs)) deallocate(dcs, stat=err)
        if (err /= 0) print *, "dcs: Deallocation request denied"
        if (allocated(drhs)) deallocate(drhs, stat=err)
        if (err /= 0) print *, "drhs: Deallocation request denied"
        if (allocated(tmp)) deallocate(tmp, stat=err)
        if (err /= 0) print *, "tmp: Deallocation request denied"
        if (allocated(b)) deallocate(b, stat=err)
        if (err /= 0) print *, "b: Deallocation request denied"
        if (allocated(residual)) deallocate(residual, stat=err)
        if (err /= 0) print *, "residual: Deallocation request denied"
        return


        999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
        call mkl_free_buffers
        stop 1
    end subroutine mso_zgmres_nonprecondition



    subroutine mso_sparse_showz(a, ia, ja, rows, cols, fmt, trans)
        use common_functions, only: prtVec
        implicit none
        complex(kind=8), dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        integer, intent(in) :: rows(2), cols(2)
        character(len=*), intent(in), optional :: fmt
        character(len=1), intent(in), optional :: trans

        integer :: i, mn(2), m, job(6), info, n
        complex(kind=8), dimension(:,:), allocatable :: adns

        m = size(ia) - 1
        n = maxval(ja)
        allocate(adns(m,n))

        job = 1
        job(4) = 2
        call mkl_zdnscsr(job,m,n,adns,m,a,ja,ia,info)

        if ( present(trans) .and. trans .eq. 'T' ) then
            do i = cols(1), cols(2)
                call prtVec( adns(rows(1):rows(2),i), fmt )
            end do
        else
            if ( present(fmt) ) then
                do i = rows(1), rows(2)
                    call prtVec( adns(i,cols(1):cols(2)), fmt )
                end do
            else
                write(*,*) adns
            end if
        end if
        if ( allocated(adns) ) deallocate(adns)
    end subroutine

    subroutine mso_sparse_showR(a, ia, ja, rows, cols, fmt, trans)
        use common_functions, only: prtVecR
        implicit none
        double precision, dimension(:), intent(in) :: a
        integer, dimension(:), intent(in) :: ia, ja
        integer, intent(in) :: rows(2), cols(2)
        character(len=*), intent(in), optional :: fmt
        character(len=1), intent(in), optional :: trans

        integer :: ii, jj, mn(2), m, job(6), info, n
        double precision, dimension(:,:), allocatable :: adns

        m = size(ia) - 1
        n = maxval(ja)
        allocate(adns(m,n))

        job = 1
        job(4) = 2
        call mkl_ddnscsr(job,m,n,adns,m,a,ja,ia,info)
        if ( .not.present(trans) .or. (trans.eq.'N' .or. trans.eq.'n') ) then
            do ii = rows(1), rows(2), 1
                do jj = cols(1), cols(2), 1
                    write(*,'('//fmt//'$)') adns(ii,jj)
                end do
                write(*,*)
            end do
        else
            if ( present(trans) .and. (trans.eq.'T' .or. trans.eq.'t') ) then
                do ii = cols(1), cols(2), 1
                    do jj = rows(1), rows(2), 1
                        write(*,'('//fmt//'$)') adns(jj,ii)
                    end do
                    write(*,*)
                end do
            else
                write(*,*) 'trans is not well defined in mso_sparse_showR.'
            end if
        end if
        if ( allocated(adns) ) deallocate(adns)
    end subroutine mso_sparse_showR




!     subroutine mso_dgmres_nonprecondition(maxiter, maxrelres, a,ia,ja,rhs, computed_solution, itercount, res,msglvl)
!         !*******************************************************************************
!         !  Content:
!         !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
!         !  msglvl: print the procesure data if msglvl>0
!         !  res(2): res(1) is the Euclidean norm of the current residual
!         !          res(2) is equal to res(1) / norm2(rhs)
!         !*******************************************************************************
!         use blas95, only: nrm2
!         implicit none
!         integer, intent(in) :: maxiter(2)
!         double precision, intent(in) :: maxrelres
!         double precision, dimension(:), intent(in) :: a
!         integer, dimension(:), intent(in) :: ia, ja
!         double precision, dimension(:), intent(in) :: rhs
!         double precision, dimension(:), intent(inout) :: computed_solution
!         integer, intent(inout) :: itercount
!         double precision, intent(inout) :: res(2)
!         integer, optional :: msglvl

!         integer :: n
!         integer :: ipar(128)
!         double precision :: dpar(128)
!         double precision, dimension(:), allocatable :: tmp
!         double precision, dimension(:), allocatable :: residual

!         integer :: rci_request, itmp

!         n = size(ia) - 1

!         allocate( residual(n))

!         if ( maxiter(2).gt.0 ) then
!             itmp = maxiter(2)
!         else
!             itmp = minval((/150,n/))
!         end if

!         allocate( tmp( (2*itmp+1)*n+itmp*(itmp+9)/2+1 ) )

!         call dfgmres_init(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
!         if (rci_request .ne. 0) go to 999
!         !---------------------------------------------------------------------------
!         ! set the desired parameters:
!         ! logical parameters:
!         ! do residual stopping test
!         ! do not request for the user defined stopping test
!         ! do the check of the norm of the next generated vector automatically
!         ! double precision parameters
!         !---------------------------------------------------------------------------
!         ipar(9) = 1
!         ipar(10) = 0
!         ipar(12) = 1
!         if ( maxiter(1).gt.0 ) then
!             ipar(5) = maxiter(1)
!         else
!             ipar(5) = minval((/150,n/))
!         end if
!         if ( maxiter(2).gt.0 ) then
!             ipar(15) = maxiter(2)
!         else
!             ipar(15) = minval((/150,n/))
!         end if
!         dpar(1) = maxrelres
!         !---------------------------------------------------------------------------
!         ! check the correctness and consistency of the newly set parameters
!         !---------------------------------------------------------------------------
!         call dfgmres_check(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
!         if (rci_request .ne. 0) go to 999
!         !---------------------------------------------------------------------------
!         ! print the info about the rci fgmres method
!         !---------------------------------------------------------------------------
!         if ( present(msglvl).and.msglvl.gt.0 ) then
!             print *, ''
!             print *,'some info about the current run of rci fgmres method:'
!             print *, ''
!             if (ipar(8) .ne. 0) then
!                 write(*,'(a,i1,a)') 'as ipar(8)=',ipar(8),', the automatic test for the maximal number of iterations will be performed'
!                 write(*,*) '(default) iteration number is ipar(5),ipar(15)', ipar(5), ipar(15)
!             else
!                write(*,'(a,i1,a)') 'as ipar(8)=',ipar(8),', the automatic test for the maximal number of iterations will be skipped'
!             end if
!             print *,'+++'
!             if (ipar(9) .ne. 0) then
!                 write(*,'(a,i1,a)') 'as ipar(9)=',ipar(9),', the automatic residual test will be performed'
!             else
!                 write(*,'(a,i1,a)') 'as ipar(9)=',ipar(9),', the automatic residual test will be skipped'
!             end if
!             print *,'+++'
!             if (ipar(10) .ne. 0) then
!                 write(*,'(a,i1,a)') 'as ipar(10)=',ipar(10),', the user-defined stopping test will be requested via rci_request=2'
!             else
!                 write(*,'(a,i1,a)') 'as ipar(10)=',ipar(10),', the user-defined stopping test will not be requested, thus, rci_request will not take the value 2'
!             end if
!             print *,'+++'
!             if (ipar(11) .ne. 0) then
!                 write(*,'(a,i1,a)') 'as ipar(11)=',ipar(11),', the preconditioned fgmres iterations will be performed, thus, the preconditioner action will be',' requested via rci_request=3'
!             else
!                 write(*,'(a,i1,a)') 'as ipar(11)=',ipar(11),', the preconditioned fgmres iterations will not be performed, thus, rci_request will not take the',' value 3'
!             end if
!             print *,'+++'
!             if (ipar(12) .ne. 0) then
!                 write(*,'(a,i1,a)') 'as ipar(12)=',ipar(12),', the automatic test for the norm of the next generated vector is not equal to zero up to rounding and',' computational errors will be performed, thus, rci_request will not take the value 4'
!             else
!                 write(*,'(a,i1,a)') 'as ipar(12)=',ipar(12),', the automatic test for the norm of the next generated vector is not equal to zero up to rounding and',' computational errors will be skipped, thus, the user-defined test will be',' requested via rci_request=4'
!             end if
!             print *,'+++'
!         end if

!         !---------------------------------------------------------------------------
!         ! compute the solution by rci (p)fgmres solver without preconditioning
!         ! reverse communication starts here
!         !---------------------------------------------------------------------------
!         1 call dfgmres(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
!         !---------------------------------------------------------------------------
!         ! if rci_request=0, then the solution was found with the required precision
!         !---------------------------------------------------------------------------
!         if (rci_request .eq. 0) go to 3
!         !---------------------------------------------------------------------------
!         ! if rci_request=1, then compute the vector a*tmp(ipar(22))
!         ! and put the result in vector tmp(ipar(23))
!         !---------------------------------------------------------------------------
!         if (rci_request .eq. 1) then
!             call mkl_dcsrgemv('n',n, a, ia, ja, tmp(ipar(22)), tmp(ipar(23)))
!             go to 1
!         !---------------------------------------------------------------------------
!         ! if rci_request=anything else, then dfgmres subroutine failed
!         ! to compute the solution vector: computed_solution(n)
!         !---------------------------------------------------------------------------
!         else
!             go to 999
!         end if
!         !---------------------------------------------------------------------------
!         ! reverse communication ends here
!         ! get the current iteration number and the fgmres solution (do not forget to
!         ! call dfgmres_get routine as computed_solution is still containing
!         ! the initial guess!)
!         !---------------------------------------------------------------------------
!         3 call dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)
!         call mkl_dcsrgemv('n', n, a, ia, ja, computed_solution, residual)
!         call daxpy(n, -1.0d0, rhs, 1, residual, 1)
!         res(1) = nrm2(residual)
!         res(2) = res(1)/nrm2(rhs)
!         !---------------------------------------------------------------------------
!         ! print the number of iterations: itercount
!         !---------------------------------------------------------------------------
!         if ( present(msglvl).and.msglvl.gt.0 ) then
!             print *, ''
!             print *,' the system has been solved'
!             print *, ''
!             print *,' number of iterations: ', itercount
!             print *,'the absolute tolerance is ', res(1)
!             print *,'the relative tolerance is ', res(2)
!         end if
!         if ( allocated(residual) ) deallocate(residual)
!         if ( allocated(tmp) ) deallocate(tmp)
!         !---------------------------------------------------------------------------
!         ! release internal intel(r) mkl memory that might be used for computations
!         ! note: it is important to call the routine below to avoid memory leaks
!         ! unless you disable intel(r) mkl memory manager
!         !---------------------------------------------------------------------------
!         call mkl_free_buffers
!         return

!         !---------------------------------------------------------------------------
!         ! release internal intel(r) mkl memory that might be used for computations
!         ! note: it is important to call the routine below to avoid memory leaks
!         ! unless you disable intel(r) mkl memory manager
!         !---------------------------------------------------------------------------
!         999 if ( rci_request.eq.-1 ) then
!             call dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)
!             call mkl_dcsrgemv('n', n, a, ia, ja, computed_solution, residual)
!             call daxpy(n, -1.0d0, rhs, 1, residual, 1)
!             res(1) = nrm2(residual)
!             res(2) = res(1)/nrm2(rhs)

!             write( *,'(a,i5)') 'this example failed as the solver has returned the error code', rci_request
!             print *, ''
!             print *,' the system has not been solved because the maximum number of iterations was reached, but the relative stopping criterion was not met'
!             print *, ''
!             print *,'current number of iterations: ', itercount
!             print *,'the absolute tolerance is ', res(1)
!             print *,'the relative tolerance is ', res(2)
!             if ( allocated(residual) ) deallocate(residual)
!             if ( allocated(tmp) ) deallocate(tmp)
!             call mkl_free_buffers
!             return
!         else
!             write( *,'(a,i5)') 'this example failed as the solver has returned the error code', rci_request
!             write(*,*) 'some errors accurred in mso_dgmres_nonprecondition!'
!             if ( allocated(residual) ) deallocate(residual)
!             if ( allocated(tmp) ) deallocate(tmp)
!             call mkl_free_buffers
!             stop 1
!         end if
!     end subroutine mso_dgmres_nonprecondition

!     subroutine mso_dgmres_precondition(maxiter, relres, a,ia,ja,rhs, computed_solution, itercount)
!         !*******************************************************************************
!         !  Content:
!         !  FGMRES ((Preconditioned) Flexible Generalized Minimal RESidual method)
!         !*******************************************************************************
!         implicit none
!         integer, intent(in) :: maxiter
!         double precision :: relres
!         double precision, dimension(:), intent(in) :: a
!         integer, dimension(:), intent(in) :: ia, ja
!         double precision, dimension(:), intent(in) :: rhs
!         double precision, dimension(:), intent(inout) :: computed_solution
!         integer, intent(in) :: itercount


!         integer :: n, m
!         integer :: ipar(128)
!         double precision :: dpar(128)
!         double precision, dimension(:), allocatable :: tmp, b, residual

!         integer :: rci_request, i
!         double precision :: dvar
!         double precision, external :: dnrm2
!         double precision, dimension(:), allocatable :: bilu0, trvec
!         integer :: ierr, err

!         n = size(ia) - 1    ! the size of the problem

!         allocate( b(n))
!         allocate( residual(n))

!         !---------------------------------------------------------------------------
!         ! Initialize variables and the right hand side through matrix-vector product
!         !---------------------------------------------------------------------------

!         call dcopy(n, rhs, 1, b, 1)

!         computed_solution = 1.d0

!         ipar(15) = min(150,n,maxiter)
!         allocate( tmp(n*(2*ipar(15)+1)+(ipar(15)*(ipar(15)+9))/2+1))

!         call dfgmres_init(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
!         if (rci_request .ne. 0) go to 999

!         ipar(8) = 1 ! if the value is not equal to 0, the dfgmres routine performs the stopping test for the maximum number of iterations: ipar(4) ≤ ipar(5)
!         ipar(15) = min(150,n,maxiter)
!         ipar(5) = ipar(15)
!         ipar(11) = 0
!         dpar(1) = relres

!         call dfgmres_check(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)
!         if (rci_request .ne. 0) go to 999

!         ipar(31) = 1
!         dpar(31) = 1.D-20
!         dpar(32) = 1.D-16
!         allocate(bilu0(size(a)), stat=err)
!         if (err /= 0) print *, "bilu0: Allocation request denied"
!         call dcsrilu0(n, a, ia, ja, bilu0, ipar, dpar, ierr)
!         allocate(trvec(n))

!         if (ierr .ne. 0) then
!             write(*,'(a,a,i1)') ' error after calculation of the',' preconditioner dcsrilu0',ierr
!               go to 999
!         end if

!         1 call dfgmres(n, computed_solution, rhs, rci_request, ipar, dpar, tmp)

!         if (rci_request .eq.0) go to 3

!         if (rci_request .eq. 1) then
!             call mkl_dcsrgemv('n',n, a, ia, ja, tmp(ipar(22)),tmp(ipar(23)))
!             go to 1
!         end if

!         if (rci_request .eq. 2) then
!             ! request to the dfgmres_get routine to put the solution into b(n) via ipar(13)
!             ipar(13) = 1
!             ! get the current fgmres solution in the vector b(n)
!             call dfgmres_get(n, computed_solution, b, rci_request, ipar,dpar, tmp, itercount)
!             ! compute the current true residual via mkl (sparse) blas routines
!             call mkl_dcsrgemv('n', n, a, ia, ja, b, residual)
!             call daxpy(n, -1.0d0, rhs, 1, residual, 1)
!             dvar = dnrm2(n, residual, 1)
!             if (dvar .lt. relres) then
!                 go to 3
!             else
!                 go to 1
!             end if
!         end if

!         if (rci_request .eq. 3) then
!             call mkl_dcsrtrsv('l','n','u',n,bilu0,ia,ja,tmp(ipar(22)),trvec)
!             call mkl_dcsrtrsv('u','n','n',n,bilu0,ia,ja,trvec,tmp(ipar(23)))
!             go to 1
!         end if

!         if (rci_request .eq. 4) then
!             if (dpar(7) .lt. 1.0d-12) then
!                 go to 3
!             else
!                 go to 1
!             end if
!             !---------------------------------------------------------------------------
!             ! if rci_request=anything else, then dfgmres subroutine failed
!             ! to compute the solution vector: computed_solution(n)
!             !---------------------------------------------------------------------------
!         else
!             go to 999
!         end if

!         3 IPAR(13) = 0
!         call dfgmres_get(n, computed_solution, rhs, rci_request, ipar, dpar, tmp, itercount)


!         if ( allocated(tmp) ) deallocate(tmp)
!         if ( allocated(bilu0) ) deallocate(bilu0)
!         if ( allocated(trvec) ) deallocate(trvec)
!         if ( allocated(b) ) deallocate(b)
!         if ( allocated(residual) ) deallocate(residual)
!         return

!         999 write( *,'(A,I5)') 'this example failed as the solver has returned the error code', rci_request
!         call mkl_free_buffers
!         stop 1
!     end subroutine mso_dgmres_precondition

    subroutine mso_dpardiso(mtype,msglvl,a,ia,ja,b,x)
    	! mtype: 2 for real(8) and symmetric positive definite
    	!        -2 for real(8) and symmetric indefinite
    	!		 11 for real(8) and nonsymmetric
    	! msglvl: 0 for no information output
    	!         1 for information output
        USE mkl_pardiso
        IMPLICIT NONE
        integer, parameter :: dp = kind(1.d0)
        integer, intent(in) :: mtype, msglvl
        integer, dimension(:), intent(in) :: ia, ja
        real(kind=dp), dimension(:), intent(in) :: a
        real(kind=dp), dimension(:), intent(inout) :: b
        real(kind=dp), dimension(:), intent(inout) :: x

        !.. internal solver memory pointer
        type(mkl_pardiso_handle), allocatable :: pt(:)
        !.. all other variables
        integer :: maxfct, mnum, phase, n, nrhs, error, nnz
        integer :: error1
        integer, allocatable :: iparm(:)
        integer :: i, idum(1)
        real(kind=dp) :: ddum(1)
        !.. Fill all arrays containing matrix data.
        n = size(ia) - 1
        nnz = size(a)
        nrhs = 1
        maxfct = 1
        mnum = 1

        allocate(pt(64))
        allocate(iparm(64))
        !..
        !.. set up pardiso control parameter
        !..
        !         allocate(iparm(64))

        do i = 1, 64
           iparm(i) = 0
        end do

        iparm(1) = 1 ! no solver default
        iparm(2) = 2 ! fill-in reordering from METIS
        iparm(4) = 0 ! no iterative-direct algorithm
        iparm(5) = 0 ! no user fill-in reducing permutation
        iparm(6) = 0 ! =0 solution on the first n components of x
        iparm(8) = 2 ! numbers of iterative refinement steps
        iparm(10) = 13 ! perturb the pivot elements with 1E-13
        iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
        iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
        iparm(14) = 0 ! Output: number of perturbed pivots
        iparm(18) = -1 ! Output: number of nonzeros in the factor LU
        iparm(19) = -1 ! output: mflops for lu factorization
        iparm(20) = 0 ! output: numbers of cg iterations

        error  = 0 ! initialize error flag
        !       msglvl = 1 ! print statistical information
        !       mtype  = -2 ! symmetric, indefinite

        !.. initialize the internal solver memory pointer. this is only
        ! necessary for the first call of the pardiso solver.

        !         allocate (pt(64))
        do i = 1, 64
          pt(i)%dummy =  0
        end do

        !.. back substitution and iterative refinement
        iparm(8) = 0 ! max numbers of iterative refinement steps
        phase = 13 ! 33 ! only solving
        call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
              idum, nrhs, iparm, msglvl, b, x, error)
        ! write(*,*) 'solve completed ... '
        if (error /= 0) then
            write(*,*) 'the following error was detected: ', error
            goto 1000
        end if

        1000 continue
        !.. termination and release of memory
        phase = -1 ! release internal memory
        call pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error1)

        !         if (allocated(iparm))   deallocate(iparm)
        deallocate(pt)
        deallocate(iparm)
        if (error1 /= 0) then
            write(*,*) 'the following error on release stage was detected: ', error1
            stop 1
        endif

        if (error /= 0) stop 1
    end subroutine mso_dpardiso

    subroutine mso_zpardiso(mtype,msglvl,a,ia,ja,b,x)
        ! mtype should be:
        !   3 for complex and structurally symmetric
        !   4 for complex and Hermitian positive definite
        !   -4 for complex and Hermitian indefinite
        !   6 for complex and symmetric
        !   13 for complex and nonsymmetric
        ! msglvl should be:
        !   0 for no showing the detail; 1 for showing the detail
        USE mkl_pardiso
        IMPLICIT NONE
        integer, intent(in) :: mtype, msglvl
        integer, dimension(:), intent(in) :: ia, ja
        complex(kind=8), dimension(:), intent(in) :: a
        complex(kind=8), dimension(:), intent(inout) :: b
        complex(kind=8), dimension(:), intent(inout) :: x

        !.. internal solver memory pointer
        type(mkl_pardiso_handle), allocatable  :: pt(:)
        !.. all other variables
        integer :: maxfct, mnum, phase, n, nrhs, error
        integer :: error1
        integer, allocatable :: iparm( : )
        integer :: i, idum(1)
        double precision :: ddum(1)
        !.. Fill all arrays containing matrix data.
        n = size(ia) - 1
        nrhs = 1
        maxfct = 1
        mnum = 1
        !..
        !.. set up pardiso control parameter
        !..
        allocate(iparm(64))

        do i = 1, 64
           iparm(i) = 0
        end do

        iparm(1) = 1 ! no solver default
        iparm(2) = 2 ! fill-in reordering from METIS
        iparm(4) = 0 ! no iterative-direct algorithm
        iparm(5) = 0 ! no user fill-in reducing permutation
        iparm(6) = 0 ! =0 solution on the first n components of x
        iparm(8) = 2 ! numbers of iterative refinement steps
        iparm(10) = 13 ! perturb the pivot elements with 1E-13
        iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
        iparm(13) = 0 ! maximum weighted matching algorithm is switched-off (default for symmetric). Try iparm(13) = 1 in case of inappropriate accuracy
        iparm(14) = 0 ! Output: number of perturbed pivots
        iparm(18) = -1 ! Output: number of nonzeros in the factor LU
        iparm(19) = -1 ! output: mflops for lu factorization
        iparm(20) = 0 ! output: numbers of cg iterations

        error  = 0 ! initialize error flag
        !       msglvl = 1 ! print statistical information
        !       mtype  = -2 ! symmetric, indefinite

        !.. initialize the internal solver memory pointer. this is only
        ! necessary for the first call of the pardiso solver.

        allocate (pt(64))
        do i = 1, 64
          pt(i)%dummy =  0
        end do

        !.. back substitution and iterative refinement
        iparm(8) = 2 ! max numbers of iterative refinement steps
        phase = 13 ! 33 ! only solving

        call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, b, x, error)
        ! write(*,*) 'solve completed ... '
        if (error /= 0) then
            write(*,*) 'the following error was detected: ', error
            goto 1000
        end if

        1000 continue
        !.. termination and release of memory
        phase = -1 ! release internal memory
        call pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, &
              idum, nrhs, iparm, msglvl, ddum, ddum, error1)

        if (allocated(iparm))   deallocate(iparm)
        if (allocated(pt))   deallocate(pt)

        if (error1 /= 0) then
            write(*,*) 'the following error on release stage was detected: ', error1
            stop 1
        endif

        if (error /= 0) stop 1
    end subroutine mso_zpardiso

end module MySparseOperator
